!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

program data_manage

   implicit none

   include 'data_manage_precision.include'
   include 'data_manage_parameter.include'

!----------------------------------------------------------------------------------------------------------------------------------|
!  ro, gp, mp, te, beta, vx, vy, vz, vv, bx, by, bz, bv
!----------------------------------------------------------------------------------------------------------------------------------|
   integer(4), parameter :: n_var_total = 11
   integer(4), parameter :: var_total_dim(n_var_total) = (/ 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1/)
   character(128), parameter :: var_total_name_list(n_var_total) =                                                                 &
      (/ 'ro', 'gp', 'mp', 'te', 'beta', 'vx', 'vy', 'vz', 'bx', 'by', 'bz' /)

   integer(4) :: ti, ci, output_method, output_time_start, output_time_end, format, output_level
   integer(4) :: ng, nx_no_bnd, nx_block, nx_total, ini_blocks, ratio, n_levels, nc, nx_lv, nc_lv, nx_tmp
   integer(4) :: nt, n_blocks, lv, p_pos, i, i_var, j_var, len_str, n_var, n_out, i_e, ip1, i_c, i_p, j, i_block
   integer(4) :: datavis_unit, datamap_unit, datavtu_unit, data_offset
   integer(4), dimension(:), allocatable :: n_blocks_all, var_dim

   logical(4) :: var_flag(n_var_total)

   real(4) :: gm, gmm1, rangex, zero_posx, cfl

   real(8) :: dx, eps, dx_lv
   real(8), allocatable :: t(:), dt(:)
   real(8), allocatable :: x(:), xh(:), x_lv(:), x_arr(:, :)

   real(output_real8_flag), dimension(:), allocatable :: ro, mx, my, mz, bx, by, bz, en, gp
   real(output_real8_flag), dimension(:), allocatable :: roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh
   real(output_real8_flag), dimension(:), allocatable :: ro_lv, gp_lv, mp_lv, te_lv, beta_lv,                               &
                                         vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv
   real(output_real8_flag), allocatable :: var(:, :, :), var_max(:), var_min(:)

   character(4) :: ti_char
   character(10) :: id
   character(128) :: datavis_path, var_name_list, datavis_name, directory
   character(256) :: header
   character(128), allocatable :: var_name(:)

   namelist /datavis_options/                                                                                                      &
      output_method,                                                                                                               &
      output_level,                                                                                                                &
      output_time_start,                                                                                                           &
      output_time_end,                                                                                                             &
      datavis_path,                                                                                                                &
      var_name_list,                                                                                                               &
      format,                                                                                                                      &
      header

!----------------------------------------------------------------------------------------------------------------------------------|

   write(*, *) 
   write(*, '(a)') '|------------------|'
   write(*, '(a)') ' Data manage starts:'
   write(*, '(a)') '|------------------|'
   write(*, *)

!----------------------------------------------------------------------------------------------------------------------------------|
!  check no reduplicate variable name
!----------------------------------------------------------------------------------------------------------------------------------|
   do j = 1, n_var_total
      i_var = 0
      do i = 1, n_var_total
         if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0) i_var = i_var + 1
      enddo
      if (i_var .gt. 1) then
         write(*, '(a, 1x, a)') ' ERROR: reduplicate variable name:', var_total_name_list(j)
         do i = 1, n_var_total
            if (index(trim(var_total_name_list(j)), trim(var_total_name_list(i))) .ne. 0)                                          &
               write(*, '(1x, a, 1x, a, 1x, a)') trim(var_total_name_list(j)), '<-', trim(var_total_name_list(i))
         enddo
         write(*, *)
         stop
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  prepare
!----------------------------------------------------------------------------------------------------------------------------------|
   datamap_unit = 10
   datavis_unit = 11
   datavtu_unit = 12

   open(datamap_unit, file = parameters_filename)
   read(datamap_unit, nml = datavis_options)
   close(datamap_unit)

   call system("mkdir -p " // trim(datavis_path))

   var_flag = .false.

   len_str = len(var_name_list)
   do i = 1, len_str
      if (var_name_list(i:i) .ge. 'A' .and. var_name_list(i:i) .le. 'Z')                                                           &
         var_name_list(i:i) = char(ichar(var_name_list(i:i)) - ichar('A') +  ichar('a'))
   enddo

   n_var = 0
   do i_var = 1, n_var_total
      if (index(var_name_list, trim(var_total_name_list(i_var))) .ne. 0) then
         var_flag(i_var) = .true.
         n_var = n_var + 1
      endif
   enddo

   write(*, '(a, i4)') ' Total number of variables is:', n_var
   write(*, '(a)') ' The following variables will be outputted for visualization:'
   write(*, *)

   allocate (var_name(n_var), var_dim(n_var), var_max(n_var), var_min(n_var))

   do i_var = 1, n_var_total
      if (var_flag(i_var)) write(*, '(1x, a8, 1x)', advance = 'no') trim(var_total_name_list(i_var))
   enddo
   write(*, *)

   var_name = ''
   var_dim = 0

   j_var = 1
   do i_var = 1, n_var_total
      if (var_flag(i_var)) then
         var_name(j_var) = var_total_name_list(i_var)
         var_dim(j_var) = var_total_dim(i_var)
         j_var = j_var + 1
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!     read parameters
!----------------------------------------------------------------------------------------------------------------------------------|
   datamap_unit = 10

   open (datamap_unit, file = trim(datamap_path) // 'params.mapout')
   read (datamap_unit, '(a)') directory
   read (datamap_unit, *) ng
   read (datamap_unit, *) nx_no_bnd
   read (datamap_unit, *) nx_block
   read (datamap_unit, *) nx_total
   read (datamap_unit, *) ini_blocks
   read (datamap_unit, *) ratio
   read (datamap_unit, *) n_levels
   read (datamap_unit, *) gm
   read (datamap_unit, *) rangex
   read (datamap_unit, *) zero_posx
   read (datamap_unit, *) cfl
   close (datamap_unit)
   
   write(*,*)
   write(*,*) 'Work space:  ', trim(directory)
   write(*,*)
   write(*,"(a15,i4)") ' x grid points=', (nx_total - ng * 2)
   write(*,*)
   
!----------------------------------------------------------------------------------------------------------------------------------|
!  read t
!----------------------------------------------------------------------------------------------------------------------------------|
   open (datamap_unit, file = trim(datamap_path) // 'nt.mapout', form = 'unformatted')
   read (datamap_unit) nt
   close (datamap_unit)

   nx_lv = (nx_total - ng * 2) * 2 ** (output_level - 1) + 2 * ng

   allocate(t(nt), dt(nt), n_blocks_all(nt),                                                                                       &
      x(nx_block), xh(nx_block), x_lv(nx_total),                                                                                   &
      ro(nx_block), mx(nx_block), my(nx_block), mz(nx_block),                                                                      &
      bx(nx_block), by(nx_block), bz(nx_block), en(nx_block),                                                                      &
      gp(nx_block),                                                                                                                &
      roh(nx_block), gph(nx_block), mph(nx_block),                                                                                 &
      teh(nx_block), betah(nx_block),                                                                                              &
      vxh(nx_block), vyh(nx_block), vzh(nx_block),                                                                                 &
      bxh(nx_block), byh(nx_block), bzh(nx_block),                                                                                 &
      ro_lv(nx_lv), gp_lv(nx_lv), mp_lv(nx_lv),                                                                                    &
      te_lv(nx_lv), beta_lv(nx_lv),                                                                                                &
      vx_lv(nx_lv), vy_lv(nx_lv), vz_lv(nx_lv),                                                                                    &
      bx_lv(nx_lv), by_lv(nx_lv), bz_lv(nx_lv))

!----------------------------------------------------------------------------------------------------------------------------------|
!  open file for x, ro, mx, en
!----------------------------------------------------------------------------------------------------------------------------------|
   gmm1 = gm - 1.0d0
   var_min = 1.0e10
   var_max = -1.0e10

   output_time_start = max(0, output_time_start) + 1
   output_time_end = min(nt - 1, output_time_end) + 1

   do ti = 1, nt
      write (ti_char, '(i4.4)') ti
      open (datamap_unit, file = trim(datamap_path) // 'info.t' // ti_char // '.mapout', form = 'unformatted')
      read (datamap_unit) t(ti)
      read (datamap_unit) dt(ti)
      close (datamap_unit)
   enddo

   do ti = output_time_start, output_time_end

      n_out = n_out + 1

      write(ti_char, '(i4.4)') ti
      datavis_name = trim (datavis_path) // 'data' // ti_char

      write(*, *) ' Time = ', t(ti)
      i_block = 0

      if (output_method .eq. 2) open(datavis_unit, file = trim(datavis_name) // '.mapvis', form = 'unformatted')

      if (output_method .eq. 3) then
         data_offset = 0
         if (format .eq. 0) then
            open(datavtu_unit, file = trim(datavis_name) // '.vtu', form = 'binary')
            write(datavtu_unit) '<?xml version="1.0"?>' // char(10)
            write(datavtu_unit) '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">' // char(10)
            write(datavtu_unit) '<UnstructuredGrid>' // char(10)
         else
            open(datavtu_unit, file = trim(datavis_name) // '.vtu')
            write(datavtu_unit, '(a)') '<?xml version="1.0"?>'
            write(datavtu_unit, '(a)') '<VTKFile type="UnstructuredGrid" version="0.1" byte_order="LittleEndian">'
            write(datavtu_unit, '(a)') '<UnstructuredGrid>'
         endif

         if (format .eq. 0) then
            open (datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.mapout', form = 'unformatted')
            read (datamap_unit) n_blocks

            do ci = 1, n_blocks
               read(datamap_unit) id
               read(datamap_unit) lv, nc, p_pos
               read(datamap_unit) dx
               read(datamap_unit) x
               read(datamap_unit) ro
               read(datamap_unit) mx
               read(datamap_unit) my
               read(datamap_unit) mz
               read(datamap_unit) bx
               read(datamap_unit) by
               read(datamap_unit) bz
               read(datamap_unit) en
               if (nc .ne. 0) cycle
               call write_vtu_header(datavtu_unit, var_name, nx_block - 2 * ng + 1, n_var, var_dim, format, data_offset)
            enddo
            close(datamap_unit)
         endif
      endif

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</UnstructuredGrid>' // char(10)
            write(datavtu_unit) '<AppendedData encoding="raw">' // char(10)
            write(datavtu_unit) '_'
         endif
      endif

!----------------------------------------------------------------------------------------------------------------------------------|
!     read n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
      open(datamap_unit, file = trim(datamap_path) // 'data.t' // ti_char // '.mapout', form = 'unformatted')
      read(datamap_unit) n_blocks
      write(datavis_unit) n_blocks

      if (output_method .eq. 1) then
         allocate(var(sum(var_dim), nx_block, n_blocks), x_arr(nx_block, n_blocks))
      endif
      if (output_method .eq. 3) then
         allocate(var(sum(var_dim), nx_block, 1), x_arr(nx_block, 1))
      endif

      do ci = 1, n_blocks
!----------------------------------------------------------------------------------------------------------------------------------|
!        readx, ro, mx, my, mz, bx, by, bz, en
!----------------------------------------------------------------------------------------------------------------------------------|
         read(datamap_unit) id
         read(datamap_unit) lv, nc, p_pos
         read(datamap_unit) dx
         read(datamap_unit) x
         read(datamap_unit) ro
         read(datamap_unit) mx
         read(datamap_unit) my
         read(datamap_unit) mz
         read(datamap_unit) bx
         read(datamap_unit) by
         read(datamap_unit) bz
         read(datamap_unit) en

         i_e = nx_block - 1

         gp = (en - (mx ** 2 + my ** 2 + mz ** 2) / ro / 2.0d0 - (bx ** 2 + by ** 2 + bz ** 2) / 2.0d0) * gmm1

         do i = 1, i_e
            xh(i) = (x(i) + x(i + 1)) / 2.0d0
         enddo

         do i = 1, i_e
            ip1 = i + 1
            roh(i) = (ro(i) + ro(ip1)) / 2.0
            gph(i) = (gp(i) + gp(ip1)) / 2.0
            bxh(i) = (bx(i) + bx(ip1)) / 2.0
            byh(i) = (by(i) + by(ip1)) / 2.0
            bzh(i) = (bz(i) + bz(ip1)) / 2.0
            vxh(i) = (mx(i) + mx(ip1)) / 2.0 / roh(i)
            vyh(i) = (my(i) + my(ip1)) / 2.0 / roh(i)
            vzh(i) = (mz(i) + mz(ip1)) / 2.0 / roh(i)
            mph(i) = (bxh(i) * bxh(i) + byh(i) * byh(i) + bzh(i) * bzh(i)) / 2.0
            teh(i) = gph(i) / roh(i)
            betah(i) = gph(i) / (mph(i) + eps)
         enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!        for base level
!----------------------------------------------------------------------------------------------------------------------------------|
         if (lv .eq. 1) then
            dx_lv = dx
            nc_lv = nc_lv + 1

            nx_tmp = (p_pos - 1) * nx_no_bnd + ng

            do i = 0, nx_no_bnd
               i_c = i + ng
               i_p = i + nx_tmp
               x_lv(i_p) = xh(i_c)
            enddo

            do i = 0, nx_no_bnd
               i_c = i + ng
               i_p = i + nx_tmp
               ro_lv(i_p) = roh(i_c)
               gp_lv(i_p) = gph(i_c)
               mp_lv(i_p) = mph(i_c)
               te_lv(i_p) = teh(i_c)
               beta_lv(i_p) = betah(i_c)
               vx_lv(i_p) = vxh(i_c)
               vy_lv(i_p) = vyh(i_c)
               vz_lv(i_p) = vzh(i_c)
               bx_lv(i_p) = bxh(i_c)
               by_lv(i_p) = byh(i_c)
               bz_lv(i_p) = bzh(i_c)
            enddo
         endif

         if (output_method .eq. 2) then
            call write_mapvis(lv, nc, dx, xh,                                                                                      &
               roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh,                                                            &
               var_flag, nx_block, ng, n_var_total, datavis_unit)
         endif

         if (nc .ne. 0) cycle

!----------------------------------------------------------------------------------------------------------------------------------|
!        calculate the maximum and minimum value
!----------------------------------------------------------------------------------------------------------------------------------|
         call extremum_calculate(var_min, var_max, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, var_flag,              &
            n_var, nx_block, ng, n_var_total)

!----------------------------------------------------------------------------------------------------------------------------------|
!        var array and coordinate array
!----------------------------------------------------------------------------------------------------------------------------------|
         if (output_method .eq. 1) then
            i_block = i_block + 1
            call data_collect(var, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, var_flag,                              &
               x_arr, xh, sum(var_dim), nx_block, n_blocks, n_var_total, i_block)
         endif

         if (output_method .eq. 3) then
            i_block = 1
            call data_collect(var, roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh, var_flag,                              &
               x_arr, xh, sum(var_dim), nx_block, i_block, n_var_total, i_block)
         endif
         if (output_method .eq. 3) then
            call write_vtu_data(datavtu_unit, var_name, real(var(:, ng : nx_block - ng, i_block)),                                 &
               real(x_arr(ng : nx_block - ng, i_block)), nx_block - 2 * ng + 1, n_var, var_dim, format)
         endif
      enddo
      close(datamap_unit)

      if (output_method .eq. 3) then
         if (format .eq. 0) then
            write(datavtu_unit) '</AppendedData>' // char(10)
            write(datavtu_unit) '</VTKFile>'
            close(datavtu_unit)
         else
            write(datavtu_unit, '(a)') '</UnstructuredGrid>'
            write(datavtu_unit, '(a)') '</VTKFile>'
            close(datavtu_unit)
         endif
      endif

      if (output_method .eq. 2) then
         call write_mapvis(output_level, nc_lv, dx_lv, x_lv, ro_lv, gp_lv, mp_lv, te_lv, beta_lv, vx_lv, vy_lv, vz_lv,             &
            bx_lv, by_lv, bz_lv, var_flag, nx_lv, ng, n_var_total, datavis_unit)
         close(datavis_unit)
      endif

      if (output_method .eq. 1) then
         call write_vtk_unstructured(trim(datavis_name) // '.vtk', var_name, real(var(:, ng : nx_block - ng, :)),                  &
            real(x_arr(ng : nx_block - ng, :)), nx_block - 2 * ng + 1, n_blocks, n_var, var_dim, format, header)
      endif

      if (output_method .eq. 3) deallocate(var, x_arr)
      if (output_method .eq. 1) deallocate(var, x_arr)
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!   write parameters
!----------------------------------------------------------------------------------------------------------------------------------|
   datavis_unit = 11
   open (datavis_unit, file = trim(datavis_path) // 'par.mapvis')
   write (datavis_unit, '(a)') 'work space: ' // trim(directory) // trim(datavis_path)
   write (datavis_unit, '(i5,t64,a)') ng, 'ng'
   write (datavis_unit, '(i5,t64,a)') nx_no_bnd, 'nx_no_bnd'
   write (datavis_unit, '(i5,t64,a)') nx_block, 'nx_block'
   write (datavis_unit, '(i5,t64,a)') nx_total, 'nx_total'
   write (datavis_unit, '(i5,t64,a)') ini_blocks, 'ini_blocks'
   write (datavis_unit, '(i5,t64,a)') ratio, 'ratio'
   write (datavis_unit, '(i5,t64,a)') n_levels, 'n_levels'
   write (datavis_unit, '(i5, t32, a)') output_real8_flag, 'output_real8_flag'
   write (datavis_unit, '(i5, t32, a)') output_level, 'output_level'
   write (datavis_unit, '(f10.5,t64,a)') gm, 'gm'
   write (datavis_unit, '(f10.5,t64,a)') rangex, 'rangex'
   write (datavis_unit, '(f10.5,t64,a)') zero_posx, 'zero_posx'
   write (datavis_unit, '(f10.5,t64,a)') cfl, 'cfl'
   write (datavis_unit, '(i5, t32, a)') nt, 'nt'
   do ti = 1, nt
      write(datavis_unit, '(es15.6, t32, a, i3, a)') t(ti), 't (ti = ', ti, ')'
   enddo
   do ti = 1, nt
      write(datavis_unit, '(es15.6, t32, a, i3, a)') dt(ti), 'dt (ti =', ti, ')'
   enddo
   write(datavis_unit, '(i5, t32, a)') n_var, 'n_var'
   write(datavis_unit, '(t4, a5, t32, a, i3, a)') (var_name(i_var), 'var_name (i_var = ', i_var, ')', i_var = 1, n_var)
   write(datavis_unit, '(es15.6, t32, a, a5, a)') (var_max(i_var), 'var_max (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
   write(datavis_unit, '(es15.6, t32, a, a5, a)') (var_min(i_var), 'var_min (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
   write(datavis_unit, '(i5, t32, a, a5, a)') (var_dim(i_var), 'var_dim (var = ', trim(var_name(i_var)), ')', i_var = 1, n_var)
   close(datavis_unit)

   deallocate(t, dt, x, xh, x_lv,                                                                                                  &
      ro, mx, my, mz, bx, by, bz, en, gp,                                                                                          &
      roh, gph, mph, teh, betah, vxh, vyh, vzh, bxh, byh, bzh,                                                                     &
      ro_lv, gp_lv, mp_lv, te_lv, beta_lv, vx_lv, vy_lv, vz_lv, bx_lv, by_lv, bz_lv,                                               &
      var_min, var_max)

   write (*,*)
   write (*,*) ' Read data complete!'
   write (*,*)
   
   stop
end program data_manage
